On the nonlinear dynamics of a piezoresistive based mass switch based on catastrophic bifurcation

This research investigates the feasibility of mass sensing in piezoresistive MEMS devices based on catastrophic bifurcation and sensitivity enhancement due to the orientation adjustment of the device with respect to the crystallographic orientation of the silicon wafer. The model studied is a cantilever microbeam at the end of which an electrostatically actuated tip mass is attached. The piezoresistive layers are bonded to the vicinity of the clamped end of the cantilever and the device is set to operate in the resonance regime by means of harmonic electrostatic excitation. The nonlinearities due to curvature, shortening and electrostatic excitation have been considered in the modelling process. It is shown that once the mass is deposited on the tip mass, the system undergoes a cyclic fold bifurcation in the frequency domain, which yields a sudden jump in the output voltage of the piezoresistive layers; this bifurcation is attributed to the nonlinearities governing the dynamics of the response. The partial differential equations of the motion are derived and discretized to give a finite degree of freedom model based on the Galerkin method, and the limit cycles are captured in the frequency domain by using the shooting method. The effect of the orientation of the device with respect to the crystallographic coordinates of the silicon and the effect of the orientation of the piezoresistive layers with respect to the microbeam length on the sensitivity of the device is also investigated. Thanks to the nonlinearity and the orientation adjustment of the device and piezoresistive layers, a twofold sensitivity enhancement due to the added mass was achieved. This achievement is due to the combined amplification of the sensitivity in the vicinity of the bifurcation point, which is attributed to the nonlinearity and maximizing the sensitivity by orientation adjustment of the anisotropic piezoresistive coefficients.


Introduction
The last decade has witnessed an ever-increasing demand for the application of microelectromechanical systems (Bhattacharyya et al. 2008;Azizi et al. 2022;Azizi et al. 2014;Pasquale and Somà 2010). Their low costs and high sensitivity have made them a very attractive option for detecting a variety of physical and chemical quantities, including very small mass, gas concentration and temperature variations (Park et al. 2012;Lin and Wang 2006;Madinei et al. 2015). Measuring very small biomedical masses such as viruses, bacteria, biomolecules, DNA, Abstract This research investigates the feasibility of mass sensing in piezoresistive MEMS devices based on catastrophic bifurcation and sensitivity enhancement due to the orientation adjustment of the device with respect to the crystallographic orientation of the silicon wafer. The model studied is a cantilever microbeam at the end of which an electrostatically actuated tip mass is attached. The piezoresistive layers are bonded to the vicinity of the clamped end of the cantilever and the device is set to operate in the resonance regime by means of harmonic electrostatic excitation. The nonlinearities due to curvature, shortening and electrostatic excitation have been considered in the modelling process. It is shown that once the mass is deposited on the tip mass, the system undergoes a cyclic fold bifurcation in the frequency domain, which yields a sudden jump in the output voltage of the piezoresistive layers; this bifurcation is attributed to the nonlinearities governing the dynamics of the response. The partial differential equations of the motion are derived and discretized to give a finite degree of freedom model based on the Galerkin method, and the limit cycles are captured in the frequency domain by using the shooting method. The or protein has always been a very challenging issue (Jafari et al. 2017;Baguet et al. 2019;Younis 2011). Various sensing mechanisms have been proposed for small mass sensing purposes; these methods include frequency shift in the resonance zone due to the added mass (Park et al. 2012;Chauhan and Ansari 2021;Joshi et al. 2019) symmetry-breaking of the vibration (Baguet et al. 2019;) and nonlinear bifurcations (Meesala et al. 2020). From the sensing type point of view, various sensing mechanisms have been applied so far which include piezoelectric Joshi et al. 2019;Azizi et al. 2017;Kumar et al. 2011), piezoresistive (Chu et al. 2018), electrostatic (Baguet et al. 2019), magnetic (Jafari et al. 2017) and nanocrystalline ZnO thin films (Bhattacharyya et al. 2008). The dynamic range enhancement of MEMS mass sensors has also been studied. Considering the geometry of the model cantilever beams with and without a tip mass have been widely studied in the literature for various applications ranging from energy harvesting (Zhang et al. 2022;Ghavami et al. 2018) to mass detection (Kumar et al. 2011;Zhao et al. 2018). From the analysis of the motion equations point of view, various approaches including numerical integration, nonlinear perturbation techniques, isogeometric analysis (IGA), (Madinei et al. 2015;Phung-Van et al. 2017;Thanh et al. 2022;Cuong-Le et al. 2022). Zhao et al. (Zhao et al. 2018) proposed a piezoelectric-based mass sensor with enhanced sensitivity due to the operation in the nonlinear bi-stable regime. Yaqoob et al. (2022) proposed a MEMS mass sensor for analyte detection using multi-mode excitation of a resonator. Wasisto et al. (2022) proposed a phase-locked loop frequency tracking system for portable piezoresistive cantilever mass sensors; their proposed model offered a mass detection in the order of ng. Toledo et al. (2019) demonstrated the potential of a piezoelectric resonator to develop a low-cost sensor for detecting microscopic masses; their model was capable of measuring mass with a sensitivity of 8.8 Hz/ng. Setiono et al. (2020) proposed an electrothermally actuated piezoresistive mass sensor to measure and monitor the changes of mass concentration of carbon nano particles in air. Their experiments were performed on two kinds of piezoresistive cantilever sensors and tipless atomic force microscopy cantilevers where the quality factor of the electrothermally actuated sensor in the in-plane operational mode was considerably higher than the out-of-plane vibration mode. Pinto et al. (2019) represented the most important metrics in the characterization of the dynamic mass sensors; they also showed that the quality factor dominates the mass sensing sensitivity as it improves from the order of some pg in atmospheric pressure to approximately 100 fg in the vacuum condition. Nayfeh et al. (2010) developed a mathematical model for a resonant gas sensor with the structure of a cantilever beam with a tip-mass exposed to electrostatic actuation; they captured the periodic orbits in the steady state by the finite difference method and applied Melnikov analysis for the detection of the homoclinic point, but did not account for the inertial and geometric nonlinearities. Stachiv et al. (2022) developed a 3D finite element model to accurately predict the resonant frequency and the corresponding mode shapes of a nano cantilever beam and the bound analyte as the added mass; they studied the impact of the size, mass, and the position of the analyte mass on the resonant frequency and the vibrational modes of the model. Despite the great efforts and studies, especially in SARS and COVID-19 virus detection (Broughton et al. 2020;Chan et al. 2019), ultraprecise bio-detection is still a challenging problem in biomedical and engineering sciences.  proposed a molecular electromechanical system (MolEMS) consisting of an aptamer probe bound to a flexible single-stranded DNA cantilever which was connected to a self-assembled stiff tetrahedral double-stranded DNS structure which enabled super sensitive detection of proteins and small molecules in biofluids. Biomass sensing, especially for disease diagnosis, is very demanding and a very challenging problem (Ihling et al. 2020) in the engineering field, especially when the measurement range is less than pg range.
As mentioned, biomass sensing in the bio-medical research field is a highly demanding topic; in this paper, we propose a super sensitive piezoresistive MEMS mass sensor/switch whose sensitivity is enhanced by taking advantage of nonlinearity and setting the orientation of the structure and the piezoresistive layer with respect to the crystallographic coordinates of the silicon such that the highest possible sensitivity is achieved. The model consists of a cantilever beam with a tip mass at the end. The system is excited in the vicinity of the catastrophic bifurcation point, once the added mass is deposited on the tip mass, due to frequency shift as a result of mass addition, the system exhibits a jump in the time response which is applied as a measure of the added mass. The nonlinear equation of the motion, which accounts for the effect of geometric and inertial nonlinearities and the nonlinear electrostatic excitation is discretized and the frequency response curves for the reduced order model are calculated using the shooting method.

Modelling
As shown in Fig. 1, the proposed mass sensor is composed of a silicon cantilever beam of length l, width b and thickness h with a tip mass of length 2l c and width b p . The coordinate system x-y-z is attached to the center of mass at the clamped end. The tip mass is excited by a DC voltage, V DC , superimposed by an AC voltage, V AC . The initial distance between the tip mass and the substrate is denoted by g. In order to compensate for the effect of the temperature change and counteract its effect on the resistance change of the piezoresistive layers, the piezoresistive layers are connected in a Wheatstone bridge configuration (Bao 2000;Fras et al. 2018;Zhao et al. 2016), as shown in Fig. 1(a). The piezoresistive layers are represented by R 1 , R 2 , R 3 , R 4 where R 1 − R 4 and R 2 − R 3 are parallel to each other and formed by diffusion or ionimplantation on the surface of the cantilever beam. The orientation of the cantilever beam with respect to the crystallographic direction <1 0 0> is denoted by and the angle of R 1 − R 4 with respect to the direction of beam is denoted by (Fig. 1b).
The equation of the motion is given as (Azizi et al. 2022;Ali and Nayfeh 2004): where w is the transverse displacement, ρ is mass density, A is cross-sectional area and E is modulus of elasticity. Here I is area moment of inertia of the cantilever beam with respect to neutral axis, J is the rotary inertia of the tip mass and t is time.
The boundary conditions are defined as (Nayfeh et al. 2010): (1) where M is the tip mass modelled as a rigid body with mass moment of inertia J M = 1 3 (M + m)l 2 c , and m is the added mass which is assumed to be uniformly distributed on the tip mass, is the permittivity of air. Introducing the following non-dimensional parameters (Azizi et al. 2016;Zamanzadeh et al. 2020), and removing the hat notation, the dimensionless governing equation reduces to Azizi et al. (2022): In Eq. (4), α 1 and α 2 are defined as: The associated boundary conditions in non-dimensional form reduce to Nayfeh et al. (2010): Equation (4) is a homogenous partial differential equation (PDE) subjected to non-homogenous boundary conditions. To solve this equation, one possibility is to discretize the equation and then numerically integrate it over time with updated shape functions in each time step which requires a huge amount of computational time. The other possibility is to apply the Galerkin discretization method based on the linear mode shapes to the extended Hamiltonian (Nayfeh et al. 2010; Phung-Van et al. 2019) which yields a nonlinear non-homogenous ordinary differential equation (ODE). The other method is to apply the Galerkin discretization to the Lagrangian (Firoozy et al. 2017) which results in a number of non-homogenous ODEs. The latter approach is adopted in this study to derive the governing reduced order model. The shape functions are assumed to be linear which satisfy the PDE with only the linear terms retained (Eq. 4) subjected to non-homogenous linear boundary conditions (Eq. 6). For this purpose, the kinetic and the potential energies are given as: Substituting the expression: into the Lagrangian ( L = K − U ), applying the Galerkin method, and then deriving the motion equations using d dt L q − L q + Ḋ q = 0 yields the following nonlinear ODE, for a single term (i.e.n = 1) subjected to the initial conditions. Hence (dropping the subscripts on q and ) where: In Eq. (10), is a dummy parameter to carry out the integral over the length of the tip mass. For the numerical solution, the phase space variables are defined (Madinei et al. 2015;Meng et al. 2022;) and the governing nonlinear ODEs are integrated over the time to get the time response. (9)

Resistivity tensor and orientation definition
Despite the isotropic and orientation-independent mechanical properties of silicon, the piezoresistive coefficient for some single crystalline semiconductors such as silicon and germanium is dependent on the orientation and accordingly anisotropic . The stress tensor ( T ) in crystalline material causes a change in the resistivity tensor ( ). The relation between the resistivity tensor and the stress tensor is given as (Bao 2000): In Eq. (11) the T i , i = 1, 2 … 6 and i , i = 1, 2 … 6 are the components of the second rank stress and resistivity tensors and are related by the piezoresistive coefficient tensor ( ) which is a tensor of fourth rank. T i and i are represented in matrix notation as: In the crystallographic coordinate system for silicon there are only three non-zero independent components for the piezoresistive coefficients which are 11 = 22 = 33 , 12 = 21 = 13 = 31 = 23 = 32 and 44 = 55 = 66 . The components of the stress and resistivity tensor in another orientation rather than the crystallographic coordinate system are determined based on the coordinate transformation of the second order tensors; however the piezoresistive coefficient tensor obeys the corresponding rule for fourthorder tensors. The components of the piezoresistive tensor for both p-Si and n-Si are given in Table 1 (Smith 1954).
As the piezoresistive layers are positioned on the surface of the structure, in most applications the plane stress condition ( T 3 = T 4 = T 4 = T 5 = 0 ) holds and accordingly the component of the resistivity tensor along the length of the piezoresistive layers (two terminals are located at either end of the terminals) reduces to: Equation (13) implies that the sensitivity along direction 1, is Δ ∕ 0 1 = 11 T 1 + 12 T 2 + 16 T 6 . Here T 1 , T 2 and T 6 are longitudinal, transverse and shear stresses in a two dimensional element in which the longitudinal direction lies along the length of the piezoresistive layer; for simplicity the corresponding sensitivity is modified as (Δ ∕ ) l = l T l + t T t + s T s where l, s, and t refer to the longitudinal, transverse and shear respectively; the components of the piezoresistive coefficients and stress tensors are computed based on determining the and angles and the tensor transformation rules corresponding to the second and fourth order tensors respectively.

Results and discussions
The geometric and mechanical properties of the studied model are given in Table 2.
Considering the coefficients of the piezoresistive tensor in the crystallographic coordinate system ( = 0 ), the orientation dependency of the piezoresistive coefficients in terms of for both p-Si and n-Si are illustrated in Fig. 2. (13) 1 = 0 1 + 11 T 1 + 12 T 2 + 16 T 6 The longitudinal, transversal and shear stress distribution for R 1 and R 4 in terms of and are shown in cylindrical coordinates in Fig. 3; here, we assume for a horizontal cross-section (a fixed ), varies in the range of 0 • to 360 • . The components of stress tensor for the pair of R 2 and R 3 can be determined by adding 90 • to the corresponding of the pair R 1 and R 4 .
As illustrated, for = 0 , the longitudinal stress l has a maximum which does not change with the variation of . As increases, l reduces while t increases. For = 90 • , l is zero while t becomes maximum. The variation of the shear stress in terms of and is given in Fig. 3c which shows that the shear stress is maximum for = 45 Each horizontal cross section in Fig. 4 shows the variation of the stress tensor components for a fixed while varies from 0 • to 360 • . As illustrated, regardless of , the maximum value for l occurs at = 0 , however the transverse and shear stresses reach their maximum values at = 90 • and = 45 • , respectively. As discussed in the previous section, the sensitivity depends on stress and resistivity tensors. Considering a p-Si and assuming a 1 MPa uniaxial tensile stress along the length of the beam, we examine the sensitivity of the R 1 and R 4 resistors to determine the most sensitive orientation in terms of and . Figures 5a and b show the corresponding sensitivity for fixed and variable and for fixed and variable , respectively. The maximum sensitivity corresponds to = 45 • for an orientation of < 110 > and = 0 • which are set as the desired orientations in the rest of this study. The frequency response curves of the cantilever beam in the vicinity of both superharmonic and primary resonances are shown in Fig. 6 for V DC = 1v and V AC = 0.1v.
Various sources of nonlinearity such as electrostatic, inertial, and geometric nonlinearity are present, however, electrostatic nonlinearity dominates the response. Electrostatic nonlinearity is in the form of quadratic nonlinearity and accordingly doubling and halving the excitation frequency mechanisms are active in the system; This ends up with super-harmonic and sub-harmonic nonlinear resonance zones on the frequency response curves of the system. Doubling of the excitation frequency, activates of the primary resonance of the system when the micro beam is excited by half of the frequency of the primary resonance. The existence of frequency doubling (in case (a) ( b) Fig. 2 Dependency of the piezoresistive coefficients in terms of , a p-Si, b n-Si of super-harmonic resonance) and frequency halving (in case of sub-harmonic resonance) in the dynamics of the system cause this nonlinear behaviour. This suggests that the amplitude of the motion should be increased to enhance the nonlinearity which accordingly produces nonlinear bifurcation points in the frequency domain; these bifurcation points are likely to exhibit super sensitivity for mass detection purposes. Figure 7 illustrates the frequency response curves for two different values of V DC and for V AC = 0.1V. By increasing V DC , the bistable solutions approach each other and accordingly the higher amplitude branch of the solution for V DC = 7V falls under the Fig. 3 The components of stress tensor distribution (for each value of on the vertical axis varies from 0 • to360 • ) a l , b t , c s one corresponding to V DC = 4V . Once the excitation frequency is increased, the lower amplitude branch of the solution loses stability through a cyclic fold bifurcation point (Azizi et al. 2022). This point is where the stable and unstable manifolds intersect, and the system undergoes a sudden jump; the aim is to operate the system close to this bifurcation to take advantage of the sudden jump which offers a super sensitive regime for mass detection. Here we excite the system in the vicinity of the bifurcation point. The reason to not to operate exactly at the bifurcation point is because of the super sensitivity to very small disturbances since the low amplitude stable  Fig. 8 might be a loss of stability and, accordingly a jump to a higher stable branch. The frequency response curves are shown in Fig. 8 for two different cases. The first case is in the absence of the added mass; whereas in the second case it is assumed that a 100 pg mass is uniformly deposited on the tip mass which shifts the frequency response curve toward the left and this means that any excitation frequency left of the bifurcation point which used to be on the lower stable manifold of the bistable solution, might remain right of the bifurcation point after the mass deposition and accordingly settle on the upper stable branch of the frequency response curve. In this section, we excite the tip mass with the frequency = 1.44 ( 13.244kHz ) corresponding to point C which is on the left of bifurcation point A ( = 1.445 ( 13.290kHz )) with the electrostatic excitation V AC = 0.1V , V DC = 4.0V . Figure 9 illustrates the time responses of both the tip of the cantilever beam and the center of the tip mass; here we assume that a uniformly added mass 10 pg is deposited on the tip mass at point C. The added mass shifts the frequency response curve to the left and accordingly, the excitation frequency approaches the bifurcation point on the modified frequency response curve. As represented here the added mass is not big enough to push the excitation frequency to the right of bifurcation Time responses before and after 10 pg added mass a x = 1 on the cantilever beam and b the center of the tip mass point A and accordingly enable the jump to the higher amplitude stable solution. Once the 10 pg mass is deposited on the tip mass, the amplitude of the motion increases slightly which is attributed to the approach of the excitation frequency to the resonance zone but not to the jump. Figure 10 gives the time responses of the cantilever tip and the center of the tip mass before and after the deposition of the 100 pg added mass.
With the deposition of the 100 pg added mass, the frequency response curve moves to the left more than the distance of the excitation frequency from the bifurcation point and accordingly, the system On the nonlinear dynamics of a piezoresistive based mass switch based on catastrophic… jumps through point D (Fig. 8) as the only possible stable response to which it eventually settles. The corresponding longitudinal normal stress l and sensitivity variations along the R 1 and R 4 resistances are illustrated in Fig. 11. We demonstrated that for the response to jump to the higher stable branch of the solution, the added mass needs to be more than the minimum required to trigger the bifurcation. It was shown that 10 pg added mass does not trigger the bifurcation in case the system is set to operate at C on the frequency domain; however, Fig. 12 shows that a jump is trigger when the system is set to initially operate at point B (shown in Fig. 8) with = 1.4445 ( 13.285kHz ) which is 4.6 Hz from the bifurcation point A.

Conclusion
In this paper, the nonlinear dynamics of a cantilever microbeam with a tip mass for mass sensing was investigated. The tip mass was subjected to electrostatic actuation as a combination of DC and AC voltages. Four piezoresistive layers in a Whetstone bridge configuration were connected to each other in the vicinity of the clamped end which undergoes the maximum normal stress. Since silicon is a non-isotropic material in terms of resistivity, we accounted for the effect of the orientation of the structure with respect to the crystallographic coordinates and its effect on the components of the resistivity tensor. The dependency of the sensitivity on the orientation of the piezoresistive layers with respect to the cantilever length was also investigated. We accounted for the geometric, inertial and electrostatic nonlinearities in the modelling. It was demonstrated that due to the nonlinearity of the response, there exists a cyclic fold bifurcation point on the frequency response curves where the stable and unstable branches of the solution on the frequency domain approach each other and they intersect at the cyclic fold bifurcation point beyond which both branches disappear and accordingly a jump occurs. It was shown that mass deposition moves the frequency response curves leftward in the frequency domain. Should the system be run to the left of the jump bifurcation, for an added mass higher than the threshold value, the operation frequency remains to the right of the bifurcation point, where system finds no stable low amplitude solutions and consequently jumps to the high amplitude stable response. We not only took the benefit of the nonlinearity to force the system to exhibit a sudden jump in the amplitude of the motion once the added mass is deposited on the tip mass, but also, we investigated the most sensitive direction of the piezoresistive layers with respect to the crystallographic orientation of the silicon. The study was performed for both p-Si and (a) ( b) Fig. 12 Time responses before and after 10 pg added mass a x = 1 on the cantilever beam and b the center of the tip mass n-Si. The joint enhancement of the sensitivity due to the nonlinearity and orientation adjustment of the piezoresistive layers, offered a twofold sensitivity enhancement which is a promising improvement in the design of mass sensors for biological applications.

Conflict of interest
The authors declare that they have no conflict of interest.
Data availability Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.